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We investigate the equilibrium and dynamical properties of the Bose-Hubbard model and the 
related particle-hole symmetric spin-1 model in the vicinity of the superfiuid to Mott insulator 
quantum phase transition. We employ the following methods: exact-diagonalization, mean field 
(Gutzwiller), cluster mean-field, and mean-field plus Gaussian fluctuations. In the first part of the 
paper we benchmark the four methods by analyzing the equilibrium problem and give numerical 
estimates for observables such as the density of double occupancies and their correlation function. In 
the second part, we study parametric ramps from the superfiuid to the Mott insulator and map out 
the crossover from the regime of fast ramps, which is dominated by local physics, to the regime of 
slow ramps with a characteristic universal power law scaling, which is dominated by long wavelength 
excitations. We calculate values of several relevant physical observables, characteristic time scales, 
and an optimal protocol needed for observing universal scaling. 



I. INTRODUCTION 

Recent experimental progress on trapping, control, and 
imaging of ultra-cold atoms has advanced to a similar 
level of accuracy as high fidelity simulation of the Bose- 
Hubbard model over a wide parameter regime. Most in- 
triguing aspect of experiments reported in Refs. [THSl is 
the ability to read out the state of the atoms with sin- 
gle site resolution. Although the static properties of the 
Bose-Hubbard model have been studied extensively using 
various numerical techniques, most significantly Quan- 
tum Monte Carlo (QMC), the same is not true regarding 
its dynamic properties. In particular, we are interested in 
the time dependent behavior of the Bose-Hubbard model 
undergoing parametric drive (i.e. tuning of the hopping 
matrix element and the onsite interaction in time). 

One aspect of parametric drive that has received a 
lot of theoretical attention is ramping a system across 
a phase transitiorPP The classical version of this prob- 
lem was originally addressed by Kibble and Zurek, who 
observed that for a thermodynamic phase transition as 
a system is driven from a disordered phase into an or- 
dered phase, different regions of the system order inde- 
pendently, thus, entrapping topological excitations (e.g. 
vortices or hedgehogs). Further, the density of topologi- 
cal excitations that is entrapped depends on how fast and 
how far into the ordered phase the system is driverP^. 
Later, this analysis was extended to dynamical crossing 
of a quantum phase transitiorP^. There, it was ob- 
served that dynamics becomes universal if one paramet- 
rically drives the system in the vicinity of a quantum 
critical point. Moreover, the signatures of this universal 
dynamics will appear as universal power law behavior of 
observables like the number of collective modes excited 
and the energy pumped into the system as a function 
of how rapidly the system is ramped across the phase 



transitiorPHH. 

The Bose Hubbard-Model and the corresponding ex- 
periments using ultracold bosons in optical lattices, are 
a natural test bed for studying universal dynamics both 
theoretically and experimentally. From the experimental 
perspective, studying bosons in lattices is attractive for 
several reasons. First, it is possible to prepare the atoms 
in very "cold" superfiuid state, thus, allowing one to con- 
centrate on studying the quantum dynamics. Second, the 
timescales available in experiments are favorable for ob- 
serving dynamics: the extrinsic timescale for atom loss is 
long compared to the intrinsic Bose-Hubbard timescales 
which are themselves long compared to the time needed 
to tune the system parameters and make observations. 
Finally, using the ultracold gas microscope^', it is possi- 
ble to see the occupation number (or at least its parity) 
of individual lattice sites, thus, obtaining a very sensitive 
probe of the dynamics. The parity probe can be thought 
of as counting defects, i.e. sites with too many or too 
few bosons as compared to the average occupation num- 
ber. From the theoretical perspective, the equilibrium 
properties of the Superfluid-Mott transition are well un- 
derstood. However, there is a lack of tools that can be 
used for studying the dynamics of interacting systems. 
The Bose-Hubbard model, thus, provides theory with a 
challenge, as well as a hope of future comparison with 
experiments. 

In this manuscript, we shall apply a number of numeri- 
cal approaches to investigate the properties of the super- 
fluid to Mott insulator transition in two dimensions, at 
commensurate filling, both in and out of equilibrium. In 
particular, we shall begin by investigating the zero tem- 
perature equilibrium properties such as density of defects 
and their correlation functions. Next, we shall look at 
the same properties in dynamic ramps across the phase 
transition starting from the zero temperature equilibrium 
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state. We shall use (1) Exact Diagonalization (ED), (2) 
Mean-Field (MF) theory, (3) Cluster-Mean-Field (CMF), 
and (4) Mean-Field theory with Gaussian fluctuations 
(MF+G) methods. The MF+G me thod is an extension 
of the normal mode analysis of RefsP^ES] to include time 
dependence of the mean field on top of which the normal 
modes are constructed. We shall apply these methods to 
the Bose-Hubbard model, as well as the closely related 
spin-1 quantum rotor model. Our goals are two-fold (a) 
to model the current generation of experiments, which 
take place on relatively short timescales and (b) to com- 
ment on what is needed to observe universal scaling in dy- 
namics experiments at correspondingly longer timescales. 
We shall take the strategy of first making the confidence 
building measure of comparing the results of these meth- 
ods with QMC in equilibrium at zero temperature. In 
building up our confidence, we obtain relatively simple 
means for computing experimentally measured quanti- 
ties including the defect densities and defect correlation 
functions. Next, we shall employ these methods to study 
the dynamics during parametric ramps to the Mott insu- 
lating phase starting with equilibrium superfluid at zero 
temperature. 

In order to study dynamics numerically, it is necessary 
to make some approximations. Indeed, our first approx- 
imation will be to simplify the Hamiltonian. We do this 
in two ways. First, we truncate the Hilbert space to three 
states per site. That is, given the average density of n 
bosons per site, we truncate the Hilbert space on each 
site to the three states corresponding to occupation by 
{hq — 1, no, 7io + l} bosons. Second, for the case of MF+G 
method, we concentrate on the case of large n , which 
reduces the Hubbard model to the quantum rotor model 
and introduces an exact particle-hole symmetry. Making 
these two approximations does not effect the universality 
class of the phase transition, and, thus, should preserve 
the universal dynamics. Moreover, as we shall show these 
approximations do not change the non-universal proper- 
ties qualitatively. Finally we remark that we will focus 
on homogeneous systems, thus, avoiding the question of 
the redistributions of atoms (and energies) in the trap. 

Our main results are as follows. For understanding 
equilibrium properties, the methods we consider all have 
strength and weaknesses. All of the methods show quali- 
tative agreement for computing non-universal properties. 
The particular tests we looked at were (1) finding the 
phase boundary between superfluid and Mott insulator, 
(2) computing the defect density, and (3) defect -defect, 
particle-hole, and particle-particle correlation functions. 
The location of the phase boundary can be computed us- 
ing MF and CMF methods. We find that both of these 
methods obtain a similar phase boundary, however the 
CMF method with large clusters provides a significant 
improvement over the MF method when compared to 
the exact Quantum Monte Carlo predictions. Likewise, 
all of the methods do a reasonable job in calculating the 
number of defects in equilibrium as a function of the tun- 
ing parameter. However, they show quantitative differ- 



ence amongst themselves. These differences shrink as 
the cluster size (for CMF method) and system size (for 
ED method) increase. Finally, we can also apply these 
methods (CMF, ED, MF+G) for calculating short range 
equilibrium correlation functions. Here, we again find 
quantitative differences, but qualitatively similar results 
between the different methods. Only the MF+G method 
can be used for computing long range correlations func- 
tions, and we use it to find the diverging length scale at 
the transition and give some estimates on the quantita- 
tive values of the <?2 function near the transition (which 
could be used for comparison with experiment). 

For understanding dynamics, we find that there are 
two regimes: fast and slow. To define the timescales, 
we first define the energy scales in the problem: J is 
the hopping matrix element and U is the on-site inter- 
action matrix element (see Section In] for the exact defi- 
nition of the models that we studyjTFor fast ramps, or 
short timescales as compared to 2ixhj J (we operate near 
the phase transition where 2ttH/J ~ 4z(2ttH/U)), short- 
range physics dominants and all methods produce qual- 
itatively similar features. In particular we see a strong 
response on timescales of (2irh/U), followed by promi- 
nent oscillations with period ~ 2(2ttH/U). This short 
timescale is associated with the on-site repulsion energy 
scale, and appears naturally from the perspective of col- 
lective modes as a strong peak in the density of states of 
both the phase and the amplitude modes that occurs near 
the lines k x ± k y = it. For slow ramps or long timescales, 
long wavelength modes become important, although the 
corresponding density of states is much smaller. We note 
that what we call the fast regime was studied within the 
MF method in RefM The only method that can capture 
non-local entanglement and correlations that we have is 
the MF+G method, which shows significant deviations 
from the other methods. Within the MF+G method we 
find a crossover to the universal power-law scaling regime 
which is not seen by MF, CMF, nor ED methods. We 
find that for ramps that start in the superfluid and end 
deep in the Mott insulator, the crossover to universal 
scaling occurs for ramps of ~ 15 * (2nh/ J c ) where J c is 
the hopping matrix element at the phase transition. For 
the experimental conditions of RefP, 87 Rb in 1360 nm 
lattice, the corresponding timescale is ~ Is. It is possi- 
ble to achieve crossover to the universal scaling regimes 
for faster ramps by starting at the QCP, but experimen- 
tally that is detrimental as it is difficult to prepare the 
system at the QCP. Finally, we show that by ending the 
ramp deep in the Mott insulator, the excitations created 
near the QCP are converted into defects which can be 
detected experimentally. 

The manuscript is organized as follows. In Section [H] 
we define the Bose-Hubbard and the quantum rotor mod- 
els and truncate the Hilbert space. We introduce the 



four methods: ED, MF, CMF, and MF+G in Section III 



We study the equilibrium properties of the Bose-Hubbard 
and quantum rotor models using these methods in Sec- 
tion IV and dynamics properties in Section |Vj Finally, 
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FIG. 1: Comparison of Mott insulator — Superfluid phase 
boundary for the (a) Bose-Hubbard model and the (b) spin- 
1 quantum rotor model, in the chemical potential-hopping 
matrix element plane, obtained using mean field (cluster size 
lxl) and cluster mean field methods (cluster size ranging from 
2x2 to 3x4) in 2D. The vertical black line in (a) indicates the 
location of the tip of the n — 1 lobe obtained from quantum 
Monte Carlo (QMC^pl. In order to capture the N = 2 Mott 
lobe, we extended the basis for the lxl and 2x2 calculations 
to include 3 particles per site. Quantitatively the extent of 
the tip of the n = 1 lobe for the Bose-Hubbard model is: 
J h /Uh = {0.043,0.047,0.052,0.054,0.061} for the {lxl MF, 
2x2 CMF, 3x3 CMF, 3x4 CMF, QMC} methods, respectively. 
For the spin-1 quantum rotor model the tip extent is: J/U = 
{0.0625,0.070,0.077,0.080} for the {lxl MF, 2x2 CMF, 3x3 
CMF, 3x4 CMF} methods, respectively. 



we draw conclusions and discuss implications of our re- 

The main 



suits for future experiments in Section VI 
body of the Paper is supplemented by three appendices 
in which we describe the details of the equilibrium and 
dynamics of the MF+G model and the connection to ex- 
periment. 



II. MODELS 



The Bose-Hubbard Hamiltonian is 
H = -J H ^2 a\a,j + ^^2 Ui ( ni _ - 1 ) - ^2 W 

(ij) i l 



where a\ and a,i are the boson creation and annihilation 
operators at site i, nt = a\ai is the boson number op- 
erator, and [ii is the site dependent chemical potential 
that describes the trap. In the homogeneous setting, this 
model supports two types of ground states: Superfluid 
and Mott insulato^. The phase diagram in the ^jl/Uh 
— Jh/Uh plane consists of a number of Mott Insulating 
lobes surrounded by a SuperfhhoP3 (see Fig. [TJ . As we 
are ultimately interested in ramping between superfluid 
and Mott phases, we will mostly consider systems with 
integer average density (however, we shall also comment 
on the effects of the trap). 

Although the properties of the Bose-Hubbard model 
are the ultimate aim of this paper, it is difficult to obtain 
any analytical results regarding dynamics for this model 
in its original form. Therefore, in order to compare vari- 
ous methods, we shall use both the Bose-Hubbard model 
and its simplified cousin the quantum-rotor like model. 
The quantum rotor model is described by the Hamilto- 



H 



(ij) 



u 



(2) 



and is closely related to the Bose-Hubbard model. In 
both cases, we shall truncate Hilbert space to three states 
per site (with the exception of the calculation of the Bose- 
Hubbard model phase diagram, where we extend the ba- 
sis in order to capture the N = 2 Mott lobe in addition 
to the N — 1 Mott lobe). That is, for the quantum rotor 
model the site basis will corresponds to the eigenstate 
of the S* operator {| — 1),;, |0)j, Similarly, for the 

Bose-Hubbard case, we use the eigenstates of the num- 
ber operator {\no — l)j, \no)i, \no + as the site-basis. 
The main ingredient that differentiates the Hubbard and 
the rotor models is that the later is symmetric between 
\hq — l)i and \uq + l)j. However, we note that in the 
limit of large Uq the two models become identical un- 
der the identification J = uqJh and U = Uh- We shall 
ignore particle-hole asymmetric terms on the basis that 
they will not effect the properties of the phase transi- 
tion. Furthermore, the asymmetric terms will be absent 
at large fillings. Regardless of the issues of the Hilbert 
space truncation and asymmetries, the low energy, long 
wavelength theory which is essential for understanding 
the scaling in the vicinity of the Superfluid-Mott Insula- 
tor transition is identical in these two models 18 . 



III. METHODS IN AND OUT OF 
EQUILIBRIUM 

The goal of this section is to establish the various ap- 
proximate methods that can be used for treating dynam- 
ics. In particular, we investigate the use of exact diago- 
nalization (ED) on small lattices, Weiss-type mean field 
theory on single sites (MF) as well as on small clusters 
(CMF), and mean field theory augmented by Gaussian 
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fluctuations (MF+G). 



A. Exact Diagonalization (ED) 

Exact diagonalization is perhaps the most straight- 
forward of the methods available for exploring dynam- 
ics. The main disadvantage of ED is the exponential 
increase of its computational complexity as system size 
is increased. 

We diagonalize both the quantum rotor model and 
the Bose Hubbard model on small lattices with peri- 
odic boundary conditions. For the case of the Bose- 
Hubbard model the total particle number commutes with 
the Hamiltonian Eq. ([I]), and therefore we work at a fixed 
total particle number. Similarly, for the case of the rotor 
model, the conserved quantity that we fix is the total S z . 
An additional consideration is the truncation of the site- 
basis. Although, for the Hubbard model, at fixed total 
particle number, the Hilbert space is already finite, we 
find that reducing the size of the Hilbert space further 
by truncating the site-basis to three elements does not 
strongly influence the observables such as the defect den- 
sity [as we demonstrate in the next Section] . On the other 
hand, for the case of the quantum rotor model, fixing the 
total S z does not fix the size of the basis. However, we 
again find that truncating the site basis to three states 
does not significantly alter the computed observables as 
compared to including more site-basis elements. 

Using ED, it is possible to look at both static and 
dynamic properties. The static properties are obtained 
by finding the eigenstate with the smallest eigenvalue of 
the Hamiltonian in the reduced Hilbert space 

-^reduced = \<t>l)(<j>l\H\(fi m }(<f> m \, (3) 

lm 

where \<j>i) are the basis vectors that make up the re- 
duced Hilbert space. Dynamics are obtained by solving 
the Schrodinger evolution equation for the wave function 
in the reduced Hilbert space 

ifr d t \ Produced (*)) = #reduced(*) I ^reduced (*)>• (4) 

As initial condition, | produced (t = 0)), we use the eigen- 
state with the lowest eigenvalue obtained for the initial 
Hamiltonian i? re duced(i = 0). 

The properties that can be studied using ED are lo- 
cal observables, as the range of correlation functions is 
limited by system size. In particular we will look at the 
defect density and next-nearest neighbor correlations (i.e. 
particle-particle, particle-hole, hole- hole), see Table [T] for 
an explicit list. It is important to note that ED cannot 
shed light on long range correlations, and therefore phase 
boundaries. 



B. Mean field (MF) and Cluster mean field (CMF) 



We use both MF and CMF methods to investigate stat- 
ics and dynamics of both models. The key advantage of 
the mean field theories over ED is that they take some 
long range correlations into account in the form of the 
order parameter. Thus CMF forms a complimentary ap- 
proach for probing all the local observables that are avail- 
able to ED, but in addition can be used to locate phase 
boundaries and probe the dynamics of the order param- 
eter. As the cluster size increases, CMF should become 
asymptotically exact, although this regime is difficult to 
reach in practice. 

To apply the CMF method, we begin by partitioning 
the system into a set of clusters, e.g. we split up the grid 
of mx x ny sites into xy clusters of size m x n. We note 
that MF can be thought of as a special case of CMF, 
with a 1 x 1 cluster. The cluster mean field Hamiltonian 
is similar to the ED Hamiltonian, but with two important 
differences. First, the particle number (S z for the quan- 
tum rotor case) inside the cluster is not conserved, so we 
must include basis vectors with different particle num- 
bers in our truncated Hilbert space. Second, instead of 
periodic boundary conditions we couple the external sites 
of the cluster to the mean field on neighboring clusters, 
thus making the MF/CMF solutions consistent. For the 
uniform case, the cluster being diagonalized must be con- 
sistent with itself, thus making the theory self-consistent. 

For the case of the quantum rotor model, the cou- 
pling to the mean field is obtained by supplementing the 
Hamiltonian of Eq. ([2| with the boundary condition 

-t ]T Sf(Sr) + h.c, (5) 

where i runs over all sites at the boundary of the clus- 
ter being diagonalized and j runs over all sites that are 
nearest neighbors of % but are external to the clusters be- 
ing diagonalized. Wc remark that for the uniform case, 
in which the cluster must be consistent with itself, we 
must be careful in counting the number of j-sites (e.g. 
corner sites still have two neighbors in 2D). Similarly, for 
the Bosc-Hubbard case, we supplement the Hamiltonian 
over the truncated Hilbert space Eq. with 

-t ]T b\(b j ) + h.c. (6) 



To study dynamics, we evolve the wave functions and 
the order parameters in a consistent way. Explicitly, con- 
sider a pair of neighboring clusters A and B. Each time 
step can be broken down into two parts. In the first part, 
we advance the wave function on each cluster in time us- 
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Model 


Spin-1 


Bose Hubbard 


defect density Pd 


(Sf? 


|(1 +Pr,-) n eodd 
| (1 — rii) n G even 


defect correlation g 2 (l) 


(s?ns? +l ) 2 


i(l + Pri)(l + Pr i+ ,) n e odd 
- Cl — Pr- 1 ) (1 — Pr , /I rin P even 


particle-particle correlation g p ~p.2(l) 


1 QZf-t 1 QZ\ QZ (A | QZ \ 

z bf [l + bi ) b i+l (l + b z+l ) 




particle-hole correlation g p -h.2(l) 


lSf(l + Sf)Sf +l (l-S? +l ) 




hole- hole correlation Qh-hfiQ-) 


I Sf{l — Sf) Sf +[ (1 — Sf +l ) 


nh,inh,i+i 



TABLE I: Explicit list of operators that correspond to observables that we compute in the Spin-1 and the Bose-Hubbard 
models. Pr is the parity operator and is +1 if the number of bosons on site i is even, and —1 if it is odd and no is the average 
filling; n Pl j = max(ajai — no, 0) is the particle number operator that counts the excess in the number of bosons on site relative 
to no, while n^, = max(no — a\a,i, 0) is the hole number operator that counts the deficit in the number of bosons. 



ing the corresponding Schrodinger equation of motion 

\^ A {t + 8)) = \^ A {t)) + ~H A {V B {t),t)\^ A {t)), (7) 
in 

|<M* + 6)) = \yj B (t)) + ~H B (*A(t),t)\tl>B(t)), (8) 
in 

where the Hamiltonian for cluster A depends on the order 
parameter in cluster B and vice versa. In the second 
part, we recompute the order parameters (corresponding 
to the operator VE") in each cluster using the new wave 
functions 

$ A (t + 5) = {,jj A {t + S)\^ A (t + 6)), (9) 
*s(< + S) = {ip B {t + *)|%n(t + 8)). (10) 



Before proceeding to the mean field plus Gaussian fluc- 
tuations (MF+G), we develop a useful parametrization 
for the uniform MF solution of the quantum rotor model 
with the site Hilbert space truncated to three state |— ), 
|0), |+). We begin by pointing out that the MF solution 
corresponds to the product wave function 



\m<t>))=l[\m<p))i=l[[ 



cos 



|0> i + ^sm(j)(| + 1), 4- j -- 1), ) 



(11) 



where 9 and <f> are variational parameters that determine 
the wave function, and cos(2</>) sin(6')/v / 2 corresponds to 
the superfluid order parameter. We have specifically not 
included any parameters to tune the weight of the |+) 
state with respect to the |— ) state as we shall use this 
parametrization exclusively at the particle-hole symmet- 
ric point (/£ = in Eq. pj)). The equations of motion 
Eqs. (p7|)- (fl0[) correspond to the cxtremum of the effective 
Lagrangian L g = (ijj(6 ', (/>)\idt — H\ip(6,(j))) with respect 
to the 9{t) and <j>{t). Carrying out the extremization pro- 
cedure, we obtain the equations of motion 



U 



- Jz cos 2 (20) cos(0), 



9 = -2Jzsin(0) sin(2</>) cos(2</>), 



(12) 
(13) 



where z is the coordination number (e.g. for 2D square 
lattice z = 4). The stationary solutions of these equa- 
tions of motion correspond to the ground state configu- 
rations. The superfluid ground state occurs for AJz > U 
while the Mott Insulator occurs in the complimentary 
regime AJz < U, and the corresponding solutions are 



U 







for U < AzJ, (14) 



9 = 0; <j> = [7/4 - Jzcos 2 (20) for U > AzJ (15) 

Although the evolution of <f> is seemingly innocuous for 
the mean-field Mott insulating ground state, we shall re- 
visit it when obtaining the Gaussian fluctuations on top 
of the Mott insulator. 
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C. Mean field with Gaussian fluctuations (MF+G) 

Having obtained the mean held, we look at Gaussi an 
excitations on top of it. Here, we follow Refs.LUiLLiD ^ 
setup the spin- wave theory of the truncated quantum ro- 
tor model. We focus on this model, as it is significantly 
more straight forward than the Bose-Hubbard model due 
to the absence of the bosonic factors. There is some addi- 



tional simplification due to working at (S z ) = (equiva- 
lent of one particle per site average filling) , which results 
in an explicit particle-hole symmetry. 

We begin by introducing the Schwinger-boson-like cre- 
ation operators {t_ ^tj ^v, 4 } which create the states 
| — 1);, |0)i, | + 1), when acting on an "empty" site. A 
more convenient representation can be obtained by ro- 
tating these operators via 







( k,i \ 










; M = e ix 


\ h,i ) 




K J 





^sin(|)e-^ 



/ cos(f)e^ 
sin(|)e^ -icos(|)e-^ 

V &~ 



1 g-# 



^sm(|)e-^ \ 
-icos(|)e-^ 



J_p-i<t> 



(16) 



/ 



so that b corresponds to creating the "mean field" state, 
while 6^ and b^ correspond to creating the two states or- 
thogonal to the mean field state. The labels a and (j> have 
been chosen to indicate amplitude and phase modes. In 
this coordinate transformation matrix M, we have left 
the phases <j) and % as free parameters, although in equi- 
librium the phases <f> and x can be set to zero. We include 
them here in order to provide enough flexibility for the 
dynamical solutions. 

The operator b\ corresponds to an amplitude excita- 
tion, while corresponds to a phase excitation, a fact 
that we can demonstrate by considering the change in 
the order parameter if we perturb the mean-field &J, in 
the bt or b\ direction. 



{b + e*b a \S+\bl + ebl) 

w -= sin(0) - %/2cos(0)Re(e) + 0{e 2 
v2 

(b + e*b^\S+\bl + ebl) 
1 



V2 



sin 



(0) -n/2cos(0/2)I me. 



(17) 



(18) 



Without applying the perturbations, the order parameter 
in equilibrium has the value ^ sin(0). Upon applying a 

perturbation in the b^ a direction, we see that the magni- 
tude of the order parameter changes, while a perturba- 
tion in the b^ induces a change of the order parameter 
phase. Thus we identifying the nature of these opera- 
tors. As we shall see, the amplitude modes will be built 
up entirely from b a and b' a 's while the phase modes from 
be/, and b Js. This separation into amplitude and phase 
sectors, that do not mix, is a feature of the quantum ro- 
tor model, and is not present in the Bose-Hubbard model 
(where the separation only occurs at small momenta, and 
at high momenta the modes become mixed). 

Having defined the basis, our next goal is to obtain 
the effective Hamiltonian up to second order in fluctua- 



tions b a and b^. Using the effective Hamiltonian, we can 
find the ground state wave function, as well as its time 
evolution. 



Before proceeding to obtain the effective Hamiltonian, 
we make the approximation that we shall avoid imple- 
menting feedback of fluctuations of the amplitude and 
phase modes back onto the mean field in both statics and 
dynamics. That is, we shall obtain the "mean field" bo, 



described by 9{t) and 4>{t) from the solution of Eq. (12) 
and ( 13 ). For the case of static Hamiltonian 9(t) and (j>{t) 



will remain fixed, while for the case of the parametrically 
tuned Hamiltonian they will evolve as a function of time. 
On top of the mean-field solution, we shall construct the 
evolution of the fields b a and b^ without any feedback to 
6{t) and 4>(t). Our approximation of not implementing 
feedback relies on the assumption that fluctuations will 
only have a small effect on the "mean field." We comment 
that implementing feedback must be done correctly using 
additional machinery such as the Popov approximation 
in order to ensures that the phase mode remains gapless 
in the superfluid phase. 



Our goal is to expand H c s to second order in b a and b^ 
operators. To do this we convert from the Schwinger bo- 
son like formalism to the Holstein-Primakoff like formal- 



ism via the replacement b ^b j — > 1 
see Ref! 15 * 16 ( We continue by obtaining each term of the 
Hamiltonian Eq. |2]), starting with the ^ term: 



U 



5>j 



u 



J2 -(1- COB (61)) 



+ cos (6)bljb a j + -(l 



-cos (8)) bl 3 b 4 



(19) 



Next, we move on to the J term, which first requires us 
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to compute the operators 



obtain the complete Hamiltonian up to quadratic order 



V2S+ = V2 + t+,jto,j) (20) 

= cos(20) sm(9)bl tj b j - cos(20) sin(6')^ J 6 QJ 
-Msin(20) (bljboj - bl,jKj) 

<20) cos(0) (bljboj+bljbj) 



cos 



cos 



~ 2i ^ t .,6a.,- - e 2l *bl, A 



ST = S 



(21) 
(22) 



In the Hamiltonian the spin raising and lowering opera- 
tors always appear in the symmetric form 

S+ S7 + Sr S+ = cos 2 (20) sin 2 (9) + S a + 5 , (23) 



where 



(cos(9) cos(20) + i sin(20)) 2 b\b\b a jb a j + h.c. 



- (3 + cos(2<9) - 2cos(40) sin 2 (6>)) b\ d b a<i + h.c. 



2 cos 2 (20) sin 2 (6) [^6^ + 6* d b a 
S. , = \-\e Ai * (1 + cos(d)) blblb^j + h.c. 



~ (1 + cos(6>)) 6* Mt+h.c. 



cos 2 (20) sin 2 (0) f&+ ^ + 6 f , ,6^ 



(24) 



(25) 



Here, we have explicitly left terms like 6q6q in the ex- 
pression for the reason that although we know the am- 
plitude of these terms (which is unity at this level of 
approximation) we do not know their phase. We will 
shortly show that the phase associated with these terms 
may be removed by a specific choice of x m transforma- 
tion Eq. (16 1, which allows us to make the substitution 



6q6q — > 1. We note that in obtaining Eqs. ( 24|25 ) we 



have dropped terms that are first order in b a and 6^ op- 
erators (and third order in 6 ,fc operators). The reason 
for doing this is that these terms involve fro.fc operators 
at finite momentum k, which vanish since we assume a 
spatially uniform mean field. On the other hand, when 
we consider dynamics of the spatially uniform (fc = 0) 
mode, the terms that are first order in b a and b$ oper- 
ators become important. Below Eq. (401, we show that 



taking these terms into account we recover the mean field 
equations of motion Eqs. (12 1 and (13 1. 



Having obtained all parts of the Hamiltonian in real 
space, we put them all together and Fourier transform to 



H, 



U 



eff 



(1 - cos(0)) 



Jz 



cos 2 (20) sin 2 (9) 



H a2 + H$ 2 . (26) 



Here, 



H, 



<t2 



2 ^ 

k 



b~-k 



Ot<T,k Pa.k 
P*,k a cr,k 



<7, — k 



(27) 



are the quadratic Hamiltonians that describe the am- 
plitude (a = a) and the phase (a — 0) modes. The 
coefficients a a ,k smd /3 CT ,fc are explicitly given by 



a a ,k = — cos(0) - 2Jzcos 2 (20) sin 2 (#) 

+ 1 [3 + cos(26») - 2 cos(40) sin 2 

t k [cos(0) cos(20) - i sin(20)] 2 e 2i<t> 
U 



0a,k 

= ^ (l+cos(0)) - Jz cos 2 (20) sin 2 (<9) 



y[l + cos(0)] 



= _^ [ l +cosW]e 2# 



(28) 
(29) 

(30) 
(31) 



where efc = 2J(cos(foc) + cos(k y )). The Hamiltonians for 
the amplitude and phase modes are diagonalized by Bo- 
goliubov transformations, as described in Appendix | A 1| 
and the resulting spectra are presented in Fig. [2j The 
general features of these spectra are as follows: (1) The 
amplitude mode is gapped everywhere in the phase dia- 
gram except at the QCP where it becomes gapless. (2) 
On the other hand, in the superfiuid phase the phase 
mode corresponds to the Goldstone mode and is there- 
fore gapless, while in the Mott phase the phase mode 
also becomes gapped. Furthermore, the (phase) ampli- 
tude mode becomes the (anti-)symmetric combination of 
the particle and hole excitations in the Mott phase, and 
since we are working in a particle-hole symmetric model 
their spectra become identical. (3) For all spectra shown 
in Fig. |2j there is the large density of states near U/2 
coming from high momentum modes, which will play an 
important role in fast dynamics. We note that the enu- 
merated features of these spectra will be preserved for 
the case of the Bose-Hubbard model, although the spec- 
tra will change quantitatively. 



1. Wave function 

For each mode a € {a, 0} and momentum fc, the cor- 
responding Hamiltonian is quadratic and can be solved 
by a Bogoliubov transform. Therefore, in terms of b a ^k 
operators the ground state for each a and k is a squeezed 
coherent state. Since the various modes and momenta 




FIG. 2: Equilibrium dispersion relations qA3J) of the phase (dashed line) and amplitude (solid line) modes at three points in 
the phase diagram (from left to right: superfluid, QCP, and Mott insulator). The position along the (J/U) axis in the phase 
diagram is measured relative to the ratio (J/U) c at the quantum critical point. In plotting these dispersions we set k y — and 
varied k x . 



are non-interacting, we can write the total wave function 
(in terms of Holstein-Primakoff operators) as a product 
of squeezed states: 



i*> = n i^.fe) 

a£{a,<f>}, k>0 



(32) 



\i>a,k) = v/lH^IV^e^U'-Ulo), (33) 

where c a ^ is the squeezing parameter, is an ad- 
ditional phase, and |0) corresponds to the vacuum of 
Holstein-Primakoff (i.e. |0) = rii b o I empty)). The 
phase £ CTj fc contributes to the overall phase of the wave 
function and is therefore unimportant. The equilibrium 
value for the wave function parameter c . k is listed in 
Appendix [A] In the same Appendix, we obtain expres- 
sions for defect density and the correlation functions in 
equilibrium, the properties of which are the subject of 
the next section. 



2. Schrodmger equation of motion 

During parametric tuning, the Hamiltonian given by 
Eq. (26 1 remains quadratic. As a result, the wave func- 



tion at all times may be written in the squeezed form of 
Eq. ( 33 ) . Therefore to understand the evolution of the 



wave function, we need to study the evolution of the pa- 
rameters c at k(t) and Ca,k(t)- The evolution of the wave 
function in each mode is governed by the Schrodinger 
equation 



ih 



d t -2A 



a,k 



IJ 



K.-k 



-k u <? 



(34) 



Here, the second term in the square brackets is a Berry 
phase term that is associated with the fact that the op- 
erator b a ^k itself transforms in time as the mean field 
evolves. The Berry phase A(t) can be obtained by looking 
at the action of the coordinate transformation M(t) from 



Eq. dTel) on the vector bj(t) = (&o,i(*)j &«,»(*)> &&i(*)): 



h(t) ->&,(< + 8) = M(t + S) ■ ti(t) (35) 
= M(t + S)-M- 1 (t)-b i (t) (36) 
= [M{t) + Sd t M(t)] ■ Ar l {t) ■ b^t), (37) 

where t(t) — (t o i ,t + i ,t_ i). Thus we obtain the Berry 
phase 



A 



[d t M(t)]-M-\t) 



b{t + 5)- bit) 
6 

/ i (y + cos(0)(£) -\Q + is\nie)4> 

\6 + ism{9)4> i(x-cos(9)(j) 



(38) 



i [X ~ </>) 



(39) 



Setting x — — cos(#)(/>, we remove the Berry phase for 



the bo operator and thus allowing us to set b\b\ — > 1 in 



Eqs. (24) and (25). Next, by the assumption that the 



mean field is uniform, we drop the off-diagonal terms 
in A that involve modes with finite k. Thus for finite 
k, the non-zero terms of the Berry phase are: A a a — 
— 2i cos(9)<j> and Aa,^ = — i(cos(0)+l)<^>. Using this Berry 
phase, we find the evolution equation for the squeezing 
parameter 

ih(dt - 2A^ a )c a . k = 2a aM c a>h + (3^ k + /3* k c 2 a k (40) 

In the above paragraph, we have focused on the dy- 
namics of finite k modes. Now we re-examine the dy- 
namics of the k = mode, described by the vector &o,i 
of Eq. (fT6| . The time evolution of the direction of bo j is 



described by the mean field equations of motion Eqs. ( 12 ) 



and (13). As we now demonstrate, the mean field equa- 



tions of motion can be expressed in terms of &o,i> b a> 
and b,p t i operators. 



Consider the equation of motion Eq. 34 for the k = 
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component of bo j. The right hand side of Eq. 34 involves 



the Hamiltonian. For the k=0 case, the most significant 
contribution to the Hamiltonian comes from the terms 
that are first order in b a ^ and b^j operators (which van- 
ished for k 7^ case by assumption of a uniform mean 
field): 

^uniform = - ^ Sm(0) (&J, )i &0,i + &0,j&a,i) 
i 

+ JzJ2 [cos 2 (20)sin(20) (tf a>i bQ,i + 

i 

+tsin(40) sin(0) {b\ M - , (41) 

where we used &o^&o,i — > 1< On the left hand side of 
Eq. |34[ we find that &o,i is coupled with b a ^ by a Berry 
phase Eq. ( 39 ) . This Berry phase is associated with the 
time evolution of the direction of &o,i) and is precisely 
canceled by the Hamiltonian on the right hand side if 
the mean field equations of motion Eqs. (12 1 and (13) 
are satisfied. 

An alternative view of time evolution can be gained by 
noting the fact that the Hamiltonian is always quadratic 
means that the Heisenberg equations of motions for 



In terms of the squeezing parameter of Eq. ( 33 1 , the de- 
fect density is 



the quadratures (b* k b a .k + b 



ubn.-l 



and 



(b a ,kb(T,-k) close on themselves (see also Ref P° | 21 1 for 
5 = 1 superfluid tuning). That is, the equations of mo- 
tion do not involve higher order terms. This fact can be 
exploited to study the evolution of any state that can 
be described by these three quadratures, including not 
only the ground state but also the thermal state. In 
fact, the evolution of the ground state using the quadra- 
tures method exactly matches the Schrodinger equation 
method Eq. (34). We summarize the quadrature method 
in Appendix B^ 



3. Defect density operator 

As we have already stated, one of the most experi- 
mentally useful observables is the defect density. In the 
truncated quantum rotor model, the defect density oper- 
ator corresponds to (Sf ) 2 , and its expectation value, up 
to quadratic order may be written in the form 

(Pd) =\ [1 - cos(0)] + cos(fl) fc b afe ) 

fc 

+ i[l + cos(0)]$>i,* 6 *.*>- ( 42 ) 



(P d )=l[l-cos(0)]+co S (0)£ r 

fc 

+ i[l + cos(0)]£ T 



\Cakr 

\c<pk\ 



w (43) 



In Appendix[XJ we provide details on how to compute the 
defect density, as well as the various correlation functions 
from Table [I] in equilibrium. The results of the equilib- 
rium calculations are presented in section|lVj The expres- 
sion for the defect density is again used in conjunction 
with the Schrodinger equation of motion of the previous 
subsection to find the evolution of the number of defects 
following a ramp of the Hamiltonian parameters. The re- 
sults of dynamics calculations are presented in section [V] 



IV. RESULTS IN EQUILIBRIUM 

All of the proposed methods (MF, MF+G, CMF, and 
ED) have advantages and disadvantages. The goal of this 
section is to establish these advantages and disadvantages 
in the simple context of a homogeneous system in equilib- 
rium. As a confidence building measure, we try out the 
various approximate methods on both the Bose-Hubbard 
model as well as the quantum rotor model. In particu- 
lar, we will compare the location of the phase boundary 
obtained using the various approximate methods to the 
one predicted by quantum Monte Carlo^. 

In addition to building our confidence, wc will find 
some simple, yet useful, expressions for experimentally 
measured quantities. In particular, we obtain some ex- 
pressions for defect density as well as the correlation 
functions listed in Table [T] Defect statistics are of par- 
ticular interest because they are accessible in current ex- 
periments^^]. In addition, in Mott insulators defects, i.e. 
sites containing too many or too few bosons, roughly cor- 
respond to quasiparticle and quasihole excitations of the 
systems. One may expect (falsely) that upon crossing 
the phase boundary from the SF to the MI, the density 
of defects must decrease sharply. However, defect density 
is largely a local property of the system, and like other 
local properties does not show any significant structure 
near the phase transition. Using CMF and ED calcula- 
tions we will show that nothing dramatic happens to the 
density of defects, nor short-range correlation functions 
listed in Tabic [I] at the transition point. However, the 
non-local physics must indeed show a diverging length- 
scale near the phase transition point, which we identify 
using MF+G method. Detection of this diverging length- 
scale could be a good way to locate the phase transition 
experimentally, therefore we derive simple closed form 
expressions for the correlation functions. 

This section is organized as follows. First, we compare 
the location of the SF-MI phase boundary obtained by 
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various methods with its true location as obtained by 
quantum Monte Carlo simulations. Next, we compute 
the defect density using the various methods. Finally, 
we move on to compute the gi correlation functions. 



A. Phase boundary 

A first test of the quality of the approximation is the lo- 
cation of the SF-MI phase boundary. Amongst methods 
that we are investigating, MF and CMF methods are able 
to find the phase boundary. The location of the phase 
boundary that can be found using MF+G method will 
coincide with the MF method, as we do not implement 
feedback of fluctuations onto the order parameter. As 
a standard, we compare the results of the various meth- 
ods to those obtained from Monte Carlo. We note that 
all MF methods will generically tend to overestimate the 
importance of the ordered state, which in this case is the 
superfluid state. The reason for this tendency is that 
mean field methods do a poor job of taking into account 
the effects of long wave-length fluctuations that tend to 
destroy the emerging ordering. This problem is especially 
exacerbated in lower dimensions, where fluctuations are 
most important. 

We start by lookingat the Bose-Hubbard model in two 
dimensions. In Fig. Ilk, we compare the QMC resultP^ 
with MF and CMF results for various cluster sizes. We 
see that the CMF method approaches the QMC result 
as clusters get larger. Next, we plot the phase boundary 
for the Spin-1 model [see Fig. [lb]. Again we see that 
the extent of the Mott phase increases as the cluster size 
increases. However, in both cases the extent of the Mott 
phase begins to saturate for 3x4 clusters, which are the 
largest clusters that we simulated. 

In conclusion, the CMF method provides reasonable, 
but not quantitatively exact, estimates of the phase 
boundary. These estimates are significantly better for 
larger clusters. On the other hand, the qualitative fea- 
tures of the phase boundary are quite reasonable. There- 
fore, we suggest that using CMF method (or ED method 
with equivalent size) to study dynamics cannot provide 
accurate answers, but it can be used for qualitative an- 
swers. Further, as the location of the phase boundary 
depends on cluster size, one should adjust J/U to com- 
pensate, especially if one wishes to study dynamics in the 
vicinity of the transition. That is, instead of using J/U 
as the tuning parameter one should use (J/U)/(J/U) C . 

B. Defect density 

Defect density, as defined in Table [I] is especially in- 
teresting as it is measured experimentally. Currently, 
experiments can only distinguish whether the number of 
bosons on a particular site is even or odd so instead of 
the full counting statistics what is available is the density 
of "defects." 
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FIG. 3: Defect density in the Spin-1 model as a function of 
J/U along the particle hole symmetric line (/i = 0.5). CMF 
and ED methods converge away from the phase transition 
as the system/cluster size is increased. On the other hand, 
close to the phase transition the convergence is much worse. 
The MF+G method seems to always overestimate the defect 
density (as compared to the converged result which should 
lie between the 3x2 CMF and the 3x3 ED). The poor per- 
formance of the MF+G method is likely due to the lack of 
self-consistency. 



In Fig. [3J we plot the defect density as a function of 
J/U along the particle-hole symmetric line for the case of 
the Spin-1 model for ED, MF, CMF, and MF+G meth- 
ods. Comparing CMF and ED results with different clus- 
ter/system sizes, we see that the the ED and CMF results 
seem to converge as the cluster/system size increases. 
The convergence is good everywhere except in the vicin- 
ity of the phase transition. The lack of convergence in 
the vicinity of the phase transition is related to the im- 
portance of long wavelength fluctuations which cannot 
be captured by small cluster/systems size methods. 

An important question is the qualitative behavior of 
the defect density near the transition point. Here, the ED 
method shows the defect density monotonically increas- 
ing with J/U without any singular features in the vicin- 
ity of the phase transition (since ED works on finite sized 
systems, this lack of singular behavior is expected). On 
the other hand MF and CMF methods also show the de- 
fect density monotonically increasing with J/U but with 
a kink at the phase transition point, while the MF+G 
shows a slight bump in the vicinity of the phase transi- 
tion. In the CMF method we see that the kink becomes 
smoother as the size of the cluster is increased. We ar- 
gue that the defect density is mostly a short wavelength 
phenomenon and therefore should be only weakly sensi- 
tive to long wavelength fluctuations. On the other hand, 
the kink/bump is associated with long wavelength fluctu- 
ations, which are given excess importance within mean 
field theory description. Thus, we expect that a kink, 
but probably not a bump, can be found in the vicinity 
of the phase transition. However, as suggested by CMF 
with larger cluster sizes, this feature could be rather sub- 
tle, and thus it is a poor way of detecting the transition 
experiment ally. 
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C. Correlation functions 

In this subsection we investigate the short range (near- 
est neighbor) correlation functions using CMF, ED, and 
MF+G methods as well as long range correlation func- 
tions using the MF+G method. Our goal is to benchmark 
these methods as well as to establish some relatively sim- 
ple results which could be of use for experimental data 
fitting. 

We begin by looking at the defect-defect correlation 
function 

ft d (i) = (P d (i + i)P d (i)) (P d (i + i))(p d (i)). (44) 

Using ED, CMF, and MF+G methods we shall compute 
the nearest-neighbor defect-defect correlation function 
ft d (l = 1) = {Pd{i + l)P d {i)) - {Pd{i + 1)){PS)) and 
compare the results. We expect that deep in the Mott 
insulator, there are no defects, and therefore f2~ d (l — 
1) — > 0. On the other hand, as we move into the super- 
fluid the bosons become weakly interacting and therefore 
the defect-defect correlations decrease (there is a sub- 
tle point that for the case of the Spin-1 model particle- 
particle and particle- hole correlations persist, even in the 
non-interacting case J/U — > oo). In the vicinity of the 
transition, the interactions are strong and the defect den- 
sity is large, thus we expect a maximum in f d ~ d (l = !)• 
Indeed, plotting f d ~ d (l = 1) [see Fig. [iji], we find that 
all methods agree qualitatively. Quantitatively, we also 
find reasonable agreement with two exceptions. First, the 
2x2 ED method does a particularly poor job, as nearest 
neighbors are doubly connected to each other - once by 
the direct bond and a second time by the periodic bound- 
ary condition. Second, the MF+G method shows a jump 
at the phase transition, the CMF method shows a kink, 
and the ED shows a smooth curve. The origins of the 
disagreement at the phase transition are similar to those 
stated for the defect density, with one exception. The 
jump in the MF+G method is somewhat artificial, as we 
have used quadratic order expression on the superfluid 
side but are forced to resort to quartic order expression 
on the Mott side as quadratic order expression becomes 
zero, see Appendix |A 3| Indeed, the jump is replaced by 
a kink if we use quartic order expression on both sides. 

To gain better understanding of the defect-defect cor- 
relation function, wc look at its constituents: 



fd - d = f p- P + f p-H 



(45) 



the particle-particle (or equivalently, at the particle- 
hole symmetric point, the hole- hole) correlation func- 
tion f~ p {l = 1) = (P p (i + l)P p (i)) - (P p (i + l)){P p (i)) 
[Fig. Hb] and the particle-hole correlation function 

fT\l = l) = (P P (i + l)Ph(i)) (P P (i + l))<iMi)> 
[Fig. [4p] . The particle-particle correlation function cap- 
tures the fact that particles tend to repel each other, 
and thus we expect it to be negative (this remains true 
even in the non-interacting case due to the hopping term 



preferring particles next to holes). On the other hand, 
we expect that the particle-hole correlation function re- 
mains positive throughout the phase diagram, capturing 
the physics of the tunneling term which tends to form 
particle-hole fluctuations on nearest neighbor sites. Us- 
ing ED and CMF methods, we find the expected trends 
in both fP' p {l = 1) and g~ h (l = 1). On the other 
hand, we find that the MF+G method predicts a region 
of J/U where /f ~ p (l = 1) > 0. We attribute this failure 
to the fact that the MF+G method only captures sin- 
gle quasiparticle physics, and thus is unable to capture 
quasiparticle-quasiparticle repulsion. 

Quantitatively, for both /f^G = 1) and f2~ h (l = 1) 
we see that the CMF and the ED method seem to con- 
verge as cluster/system size is increased. As in the case 
of defect-defect correlations, 2x2 ED does a particu- 
larly poor job due to the nearest neighbors being doubly 
connected. Near the transition, the particle-hole corre- 
lations are largely responsible defect-defect correlations, 
thus fV h {l = 1) [Fig. fy] is very similar to f^ d (l = 1) 



Fig. [4p] (up to a factor of two due to the definitions) 
and all three methods ED, CMF and MF+G give rea- 
sonable results. Finally, wc mention that we present 
the normalized versions of these correlation function, e.g. 
gt d (l = 1) = ft d (l = l)/(^) 2 + 1» ^ Appendixg 

Having verified that the MF+G method provides re- 
sults similar to ED and CMF for correlations on nearest 
neighbor sites, we move on to the question of how corre- 
lations decay at large distances. This question we attack 
using only the MF+G method, as the other two methods 
are not easily extended to the long distances. Specifically, 
we look at the decay of defect-defect correlations f d ~ d (l) 
in the vicinity of the transition, on both the superfluid 
and the Mott insulating side, Fig. [5] To compute the 
correlator we have used both the second order (in fluctu- 
ations) term given by Eq. A18 and the fourth order term 
given by Eq. \A20\ 

On the superfluid side, both the second and the fourth 
order terms are non-zero. Naively, since the second or- 
der term is non-zero, we may consider stopping at this or- 
der. However, because the second order term comes from 
fluctuations of amplitude modes only, while the fourth 
order term has contributions from both amplitude and 
phase modes, it is important to include the fourth order 
term. The reason why phase modes become important is 
that their gap closes at small momenta, while amplitude 
modes are always gapped. As a result, the contributions 
from amplitude modes show exponential decay with dis- 
tance, as can be seen by plotting the second order term: 

f d ~ d(2) (l) ~ exp(-Z/£) (dotted traces in Fig.|]). More- 
over, as we approach the transition, the length scale £ 
diverges (dotted traces become flatter in Fig. [5]). On 
the other hand, contributions from phase modes show 
power law decay with distance, which becomes dominant 
when plotting the sum of second and fourth order terms: 

ft d ( l ) ~ ft d(2 \ l ) + ft d{4) ( 1 ) ~ l ' 2 (solid traces in 
Fig. [5]) . We note that if we were to proceeding to higher 
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FIG. 4: (a) Correlation functions (a) / 2 d " d (Z = 1) = {P d (i + l)P d (i)) - {P d (i + l))(P d (i)), (b) g~ P {l = 1) = {Pp(I + l)-Pp(i)} - 
(P p (i + l))(P p (i)) and (c) / 2 p ~ h (Z = 1) = (P p (i + l)P h (i)) - (P p (i + l))(P h (i)) as a function of J/U, computed using Mean Field 
with Gaussian fluctuations, Exact Diagonalization, and Cluster Mean Field. For the curve labeled MF+G, we used quadratic 
order on the Superfluid side and quartic on the Mott side; for MF+G' we used quartic order on both sides. The MF+G curve 
shows a jump indicating the importance of quartic fluctuations for computing nearest neighbor correlations. Overall, we see that 
all methods show qualitative similarities, and tend to converge as cluster size/system size is increased. The biggest differences 
occur in the vicinity of the phase transition where MF+G, MF+G' and CMF methods show a kink and exact diagonalization 
does not. 



orders, we would expect to find more negative power- 
laws, and thus we stop at the fourth order. 

The Mott insulating side is more straightforward for 
two reasons: (1) the second order term vanishes, and 
(2) both phase and amplitude modes become everywhere 
gapped, thus f^ d (l) ~ exp(™//£), Fig. [5] Again, the 
length scale xi can be seen diverging as one approaches 
the transitions as the traces in Fig. [5] become flatter. 

We can extract the correlation length scale £, that di- 
verges at the phase transition, by fitting f d ~ d (2)2(l) on 
the superfluid side and / d_d (4)2(0 on the Mott side with 
the form ~ exp(— Z/£), Fig. [6] We find that the scaling of 
£ ~ l r l _!y i s consistent with MFT result v = 0.5, where 
r = J/U — (J/U) c is the detuning away from the QCP. 
Using these notions, we obtain an approximate expres- 
sion for ~ d by fitting Fig. [HJ which is valid near the 
phase transition, but is also reasonably accurate through 
the phase diagram 
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V. DYNAMICS OF PARAMETRIC ALLY 
RAMPED SYSTEMS 

The goal of this section is to understand the dynamics 
of the homogeneous spin-1 quantum rotor model under- 
going a parametric ramp. That is, we consider the dy- 
namics as we ramp the parameters J and U in time in the 
vicinity of the quantum critical point (QCP) separating 
the superfluid and the Mott insulator. We note that the 
fj, term in Eq. ([2| commutes with the other two terms, 
and therefore the tuning of the chemical potential has 
no effect on dynamics. We investigate dynamics across 
all timescales, from very fast ramps (timescale 1/U) to 
very slow ramps (timescale 1/ J) , with the goal of seeing 



the non-universal dynamics for fast ramps cross-over to 
critical scaling dynamics for slow ramps. We focus on 
two types of ramps (1) ramps that start in the superfluid 
phase and stop in the Mott insulating phase; and (2) 
ramps that start on the QCP and stop in the Mott insu- 
lating phase. The primary purpose of investigating the 
crossover into the universal scaling is to gain quantitative 
understanding of whether this regime can be observed ex- 
perimentally. That is we want to obtain (1) a reasonable 
estimate for how slow one needs to ramp in order to be in 
the universal regime; and (2) an understanding of what 
observable to measure experimentally and whether it is 
detectable. Indeed, we will show that the universal scal- 
ing regime sets in for ramps of timescale ~ 10/ J, and 
the density of defects (i.e. density of doubly occupied and 
empty sites) can be used as an experimental observable if 
the ramp goes sufficiently deeply into the Mott Insulator. 

We remark that ramps we consider always go from the 
ordered (i.e. superfluid) to the disordered (i.e. Mott in- 
sulating phase) state, and thus can be thought of as oppo- 
site of the usual Kibble-Zurek process that describes the 
appearance of excitations and long range ordeJ^ as the 
system is ramped from the disordered phase into the or- 
dered phase (i.e. Mott insulator — > SF). Ramps toward 
the disordered phase are technically easier to describe 
because they do not require modeling spontaneous sym- 
metry breaking via Spinoidal decomposition, which falls 
outside the realm of mean field theories (i.e. MF, CMF, 
and MF+G). 

Explicitly, in a ramp towards the ordered phase collec- 
tive modes of the system that are associated with order- 
ing become unstable. In the language of susceptibilities, 
the susceptibility Xq(t) associated with ordering 



dt'xS-t')^') 



(48) 



acquires a pole with positive imaginary frequency. Al- 
though this type of poles, with a positive imaginary fre- 
quency, appear in mean field theories there are no flue- 
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Ramp Type 


Functional Form 


SF -> Shallow MI 


V(t) = 11 + (16 - ll)t/*max 


OT7 1 , | \ TV o 

bl 1 — > Deep Ml 


K(t) = 10 + (25 - 10)f/*max 


QCP -> Deep MI (MF, MF+G) 
QCP -> Deep MI (CMF, ED) 


V(t) 
V(t) 


= 12.0315 + (25 - 12.0315)t/t max 
= 11.5949 + (25 - 11.5949)t/*max 


SF -> SF 


K(t) = io + (ii-io)VW 



TABLE II: Functional forms of the optic lattice intensity as a function of time used for the various types of ramps. Note that 
in the ramps originating at the QCP, the functional form depends on the location of the QCP. For the case of ED, which lacks 
a QCP, the location of the QCP from CMF was used. 
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FIG. 5: Decay of the correlation function f£~ d (l) = (Pd(i + 
l)Pd{i)) — {Pd(i + l)){Pd(i)} as a function of separation dis- 
tance I plotted on a semi-log scale for various values of J/U 
in the vicinity of the transition (J/U) c = 0.0625. On the su- 
perfluid side, we plot all contributions to <?2(0 — 1 including 
both second and fourth order terms with solid lines, and the 
part from second order term only with dotted lines. On the 
Mott side, the second order term vanishes, and we only plot 
the contributions from the fourth order terms. See text for 
details. 



tuations to seed their growth. On the other hand, in the 
ramp towards the disordered state, both the order pa- 
rameter dynamics and the dynamics of the other modes 
of the system are well defined within mean field theo- 
ries (i.e. MF, CMF, and MF+G). The order parameter 
starts out finite and its dynamics correspond to the de- 
cay of its amplitude. Further, all of the normal modes of 
the system remain stable, and thus well described by the 
quadratic theory. 



FIG. 6: Diverging coherence length as a function of distance 
to the phase transition on a log-log plot. The scaling is con- 
sistent with MFT expectations of v — 0.5. 



Before proceeding with the numerical calculations, 
we argue that the dynamics in the vicinity of the 
superfluid-Mott insulator QCP has two regimes: fast 
(non-universal) and slow (universal) regime. The fast 
regime, which has been the focus of the recent experi- 
ments^l, is associated with timescales comparable to the 
inverse of the bandwidth of the phase and amplitude 
modes, that is t ramp ~ 2/U ~ 1/2 J z. In this regime 
dynamics is largely local and our CMF and ED approxi- 
mations work well. On the other hand, the slow regime is 
associated with universal long-wavelength physics in the 
vicinity of the QCP, and timescales comparable to the 
inverse of the energy scale over which the dispersion is 
linear, i.e. t ram p ~ 1/J. Ramping through the Superfluid 
— ¥ Mott transition results in the creation of a number of 
excitations (phase and amplitude modes) in the system 
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FIG. 7: Defect density (immediately following the end of the ramp) as a function of the ramp time. The calculations were 
performed for four different ramp types (as indicated by row headings), using the four different methods (as indicated by column 
headings). The parameters used in calculations correspond to Rb atoms in an optical lattice, see Appendix [D] for details. Note 
that the value of the mean field (for the MF, CMF, and MF+G methods) for the QCP — > deep MI ramp is identically zero 
throughout the ramp. As a result, for the MF method there are no defects, while for the CMF and MF+G methods defects 
appear due to fluctuations. 



having a density n cx . n cx is controlled by the ramp rate, 
and is believed to exhibit a universal power law in the 
slow regime 



(49) 



where r = 1 /t ram p is the ramp rate, d is the dimensional- 
ity of the system, v is the coherence length critical expo- 
nent, and z is the dynamical critical exponent!^. Plugging 
in the values of the critical exponents for the Superfluid- 
Mott transition at integer filling of the Bose-Hubbard 
model (i.e. at the particle-hole symmetric point) d = 2, 
z = 1, and v = 1/2 we obtain rt cx ~ r 2 / 3 (at least using 
the mean field exponents - which should correspond to 
our mean-field description of the QCP). These arguments 
will be made more precise as we analyze the numerical 
data. 

To study the creation of excitations and defects upon 
ramping through a phase transition, we must specify the 
protocol for parametrically tuning J and U. One impor- 
tant consideration is that only the density of defects (i.e. 
((S z ) 2 )), but not the density of single quasi-particle ex- 
citation (i.e. phase and amplitude modes of Eq. (Al)), is 



experimentally accessible. Therefore, we will look at both 
the density of defects and the density of quasi-particle 
excitations. We note that the density of defects is not a 
constant of motion, and will fluctuate following the end 
of the ramp. Thus, for concreteness, we choose to fo- 
cus on the density of defects immediately following the 
end of the ramp. One way to make a direct connection 
between defects and excitations is to ramp the system 
deep into the Mott regime, where defects indeed corre- 
spond to excitations. Following this consideration, we 
shall consider both (1) shallow ramps that end in the 
Mott phase but close to the superfluid-Mott transition, 
and (2) deep ramps that end deep in the Mott phase 
(close to J = 0). Another important consideration, is 
the initial point for the start of the ramp. As we shall 
show, starting right at the phase transition seems to be 
advantageous for observing universal scaling behavior, as 
there are fewer fluctuations and the scaling sets in for 
shortest ramp times (~ 1/^7) • However, experimentally 
preparing the system at the phase transition is a difficult 
task, due to the long equilibration times. A summary of 
ramp profiles that we investigate is provided in Table \U\ 
and the connection between the optical lattice intensity 
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0.75 



Quasiparticle Energy/U 



FIG. 8: Density of states of the phase (upper panel) and 
amplitude (lower panel) modes, as a function of the energy 
(measured in units of the interaction parameter U) for the 
Spin-1 model (from the MF+G method). The density of 
state curves are plotted for three values of (J/U)/(J/U) C 
as indicated by the number next to the curve (consecu- 
tive curves were displaced vertically for clarity). The val- 
ues of (J/U)/(J/U) C correspond to the superfluid phase 
({J/U)/(J/U) C = 1.1), the QCP ({J/U)/{J/U) C = 1.0), and 
the Mott phase ((J/U)/ (J/U) c = 0.9). The density of state 
curves have several interesting features. First, is the appear- 
ance of a gap in the amplitude mode on both sides of the QCP. 
Second, is the appearance of a gap in the phase mode on the 
Mott side but not the superfluid side of the QCP. Finally, 
the presence, throughout the phase diagram, of a logarithmic 
singularity in the density of states, associated with momenta 
near k x ±k y — ±7r, in the vicinity of ~ U/2 for both amplitude 
and phase modes. We note that this singularity is associated 
with the fast periodic oscillations seen in dynamics [Fig. [jj . 



and the Hubbard parameters is provided in Appendix |P| 

We begin by looking at the number of defects that 
are created by fast ramps with ramp times from zero 
to 50ms ~ 2Tih/J c . As we sweep through a range of 
values of J and U, for defmitiveness we shall alway com- 
pare time scales to the value of J and U at the critical 
point: 2ith/J c = 62 ms and 2irh/U c — 3.88 ms. Naively, 
we would expect that in order to create or remove a de- 
fect we must move an atom from one site to another, 
which is associated with the tunneling time 2nh/J c . In 
Fig. [7j we plot the number of defects created as a func- 
tion of the ramp time, for various ramp profiles, calcu- 
lated using the four different methods (MF, CMF, ED, 
and MF+G). Surprisingly, we see quite a lot of structure 
in the number of defects created even for ramp times 
as short as 4 ms ~ 2irh/U c . The emergence of defects for 
such short timescales can be understood by considering a 
two site quantum rotor model with total S z = (i.e. pop- 
ulated by two bosons). The spectrum of this model, for 
small J, consists of a lower branch (composed of mostly 
|0, 0) state, with a small admixture of |1, —1) and | — 1, 1) 
states) and a pair of upper branches (composed mostly of 
|1, — 1) and | — 1, 1) states with a small admixture of |0, 0) 
state). The gap between the upper and lower branches 
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FIG. 9: Quasiparticle energy at the bottom of the phase and 
amplitude bands (k = {0, 0}) and at the peak of the density of 
states (k — {0, n}, see Fig. [8| as a function of J/U across the 
Mott insulator (LHS) to superfluid (RHS) phase transition 
for the Spin-1 model (from the MF+G method). 



is set by the on-site repulsion U. Hence, the important 
timescale for characterizing adiabaticity of the process is 
set by 2ttH/U c , and the dynamics will be non-adiabatic 
for ramps that are fast compared to this timescale. In 
the many-site system, there is no such gap, but there is 
a large density of states in the vicinity of U/2 scale [see 
Fig. . The energy of this maximum, as a function of 
(J/Uj/(J/U)c for both the amplitude and phase mode is 
plotted in Fig. |9j confirming that in the vicinity of the 
quantum critical point the maximum in the density of 
states is always near U/2. This maximum results in the 
appearance of oscillatory features at the corresponding 
timescales. 

One may wonder whether other features in the density 
of states are reflected in the dynamics. A particularly 
interesting feature is the gap that appears in the am- 
plitude mode on either side of the phase transition. In 
particular, on the superfluid side, the amplitude mode is 
associated with the Higgs phenomenon, and the gap is 
related to the Higgs mass. We remark that we find no 
direct signatures of the Higgs mass (using all of our meth- 
ods) in ramps that go all the way across the phase tran- 
sition. However, the Higgs may be excited by performing 
small ramps or quenches on the superfluid side. Alterna- 
tively, the Higgs strongly couples to lattice modulations, 
and can be studied via l attice modulation spectroscopy 
as suggested in Ref! 16 * 22 ( Experimental and theoretical 
investigations of the Higgs will be the subject of an up- 
coming manuscript!^. 

Thus far, our description of fast dynamics has not ap- 
pealed to the presence of a phase transition. Indeed, 
because we are probing features at very short timescales 
and thus high energies, the properties of the phase tran- 
sition, and in fact its very presence, are hard to see. As 
an example, we consider a set of ramps that stay com- 
pletely on the superfluid side, see last row of Fig. [7j We 
see that the timescales for the number of defects created 
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qualitatively show features similar to ramps that cross 
the phase transition. 

Having understood fast timescales, we move on to the 
slow timescales. Here, we resort to using the MF+G 
method which is able to capture the long wavelength fluc- 
tuations. Since the MF+G method involves excitations 
of single particle modes, we can study both the number 
of excitation, as well as the number of defects created 
in the phase transition. In Fig. [10] we plot the num- 
ber of single particle excitations created as a function of 
the ramp rate (inverse of ramp time), for various proto- 
cols. We see that in all cases, for slow ramps, the am- 
plitude mode displays the expected 2/3 power law. Fur- 
thermore, with the exception of the protocol in which we 
start at the QCP, the majority of excitations are carried 
by the amplitude mode. However, timescales at which 
the scaling becomes easily discernible depend on the pro- 
tocol. For the two protocols starting in the superfluid, 
the excitation number exhibits some oscillations (with 
the period ~ 2irh/ J c at criticality) before the power law 
can be observed. As a result, the power law only be- 
comes apparent for relatively slow ramp rates, starting 
with ~ 0.01ms -1 ~ ^(J c /2ttH). On the other hand for 
ramps that start at the QCP, the amplitude and phase 
modes are identical (amplitude [phase] mode corresponds 
to the symmetric [anti-symmetric] combination of parti- 
cle and hole excitations of the Mott insulator), and the 
power law scaling sets in for relatively faster ramp rates 
~ 0.1ms -1 ~ 6(J c /2irh). Both of these timescales are 
comparable to what is possible in present experiments. 

To understand the feasibility of observing power-law 
scaling in experiments, we investigate the defect den- 
sity. In Fig. [TT] we plot the components of the defect 
density due to the mean-field, the amplitude modes and 
the phase modes (i.e. corresponding to the various terms 
of Eq. (A7)). We begin by looking at SF — > shallow 
MI ramps. For these ramps the number of defects cre- 
ated quickly saturates, the saturation values correspond- 
ing to the background defect density in the Mott insu- 

In order to observe scaling 



lator [panel (a) of Fig. 11 



within this protocol, it is necessary to subtract of this 
background defect density. An alternative approach is 
to perform ramps that go deep into the Mott insulating 
region and thereby convert excitations to defects. The 
results of this approach are illustrated in panel (b) of 
Fig. |11[ which shows nice power law scaling in the de- 
fect density down to extremely long ramp times (~ 3 s). 
However, the 2/3 power law scaling associated with the 
amplitude mode is partially obscured by the defects as- 
sociated with the decaying mean field, and only becomes 
apparent for timescales Is ~ 15 * (2irh/J c ) and defect 
densities of 0.1%. We remark that all issues: (1) back- 
ground defect density, (2) oscillations, (3) obscuring by 
defects associated with the mean field and (4) low defect 
density in the scaling regime can be avoided by starting 
the ramp at the QCP and going deep into the Mott insu- 
lator [panel (c) of Fig. [TT] . Using this protocol, the power 
law scaling sets in at <~ i(27r/j/J c ) and defect density of 



1%. 

Due to convenience for the MF+G method, we have 
attacked the problem of dynamics in the setting of the 
quantum rotor model. However, we expect that all the 
qualitative features will transfer directly to the Bose- 
Hubbard model, as was the case for equilibrium prop- 
erties. Thus in parametric ramps of the Bose-Hubbard 
model we expect to find the fast and slow regimes. In the 
fast regime, we again expect to see response and oscilla- 
tions on timescale of 2 * (2tvH/Uh,c)- While in the slow 
regime we expect to see a crossover to power law scaling 
at timescales of 15 * (2-kK/ Jh,c) and defect densities of 
0.05-0.5%. 



VI. DISCUSSION 

In this Paper, we have investigated the static and dy- 
namic properties of the Bose-Hubbard model, and the 
related Spin-1 quantum rotor model using four types of 
methods: (1) Mean Field, (2) Exact Diagonalization, (3) 
cluster Mean Field, and (4) Mean Field + Gaussian fluc- 
tuations. These methods were chosen as they are suitable 
for studying both statics and dynamics. We show that 
these methods are reasonably consistent with each other 
and quantum Monte Carlo calculations where available. 
The main conclusions are as follows. 

For the static case, we can compare these methods 
to Quantum Monte Carlo simulations. We find that al- 
though the qualitative features and trends displayed by 
these methods are very good, it is difficult to achieve 
high quantitative accuracy as it requires the use of large 
system/cluster sizes. In particular, we have applied the 
methods to the calculation of several observables. These 
calculations are much simpler than the corresponding 
quantum Monte Carlo calculations and thus could be of 
interest to experimentalist for data fitting. We find the 
following properties. (1) The location of the phase transi- 
tions is accessible to the various mean field methods (MF, 
CMF, and MF+G) but not to the exact diagonalization. 
The single site mean field makes an errors of 30%, by 
going to a 3 x 4 cluster the error is reduced to 11%. (2) 
The defect density in the vicinity of the phase transition 
varies smoothly within exact diagonalization, but shows 
a kink within MF and CMF methods. As the defect den- 
sity is associated with a local observable, we believe that 
any singular behavior at the phase transition should be 
strongly muted. Indeed, we see the weakening of the kink 
within CMF theory as the size of the cluster is increased. 
Further, the defect density obtained by the CMF method 
approaches the one obtained by the ED method as the 
size of the cluster (for CMF) and system (for ED) is in- 
creased. (3) Particle-particle correlation functions show 
signatures of repulsion on nearest-neighbor sites while 
particle-hole correlation functions show attraction. How- 
ever, these correlation function strongly depend on the 
method being used. (4) Finally, at large distances, the 
particle-particle correlation function shows a diverging 
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FIG. 10: Excitation density as a function of ramp rate (l/t Taulp ) on a Log-Log plot, (a) Starting in SF ending in shallow 
MI. (b) Starting in SF ending in deep MI. (c) Starting at the critical point and ending deep in the MI (here, the MF order 
parameter is zero throughout). Note that starting at the QCP results in the observation of scaling at for shorter ramp times. 
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FIG. 11: Defect density (immediately following the end of the ramp) as a function of the ramp rate (l/£ r am P ) on a Log-Log 
plot, (a) Starting in SF ending in shallow MI. (b) Starting in SF ending in deep MI. (c) Starting at the critical point and 
ending deep in the MI (here, the MF order parameter is zero throughout). Note that for trace (a), there is a large number of 
background defects in equilibrium as the ramp only goes into shallow Mott regime, resulting in the truncation of the power 
law. On the other hand, ending the ramps in the deep Mott regime results in essentially the same trace as the number of 
excitations. 



correlation length, described quantitatively by Eq. (47) 
in the vicinity of the QCP. 

For dynamics, we find two regimes: a fast regime dom- 
inated by nearest neighbor physics which is addressed 
quite well by CMF/ED methods, and a slow regime dom- 
inated by long wavelength excitations, which is addressed 
by MF+G method. The methods we use, unlike QMC, 
are well suited for treating dynamics. Following the in- 
tuition gained by comparing them in equilibrium with 
each other and with QMC, we suspect that the methods 
can be used as a guide for understanding experiments. 
Specifically, for properties like the defect density, we ex- 
pect the results to be accurate to within a factor of two. 

We find that the fast regime is dominated by the 
peak in the density of states of single quasiparticle 
modes in the vicinity of <7/2, which correspond to short- 
wavelength excitations, and results in timescales of ~ 
2 * (2irh/U c ). As a result, the short time behavior can 
probably be accurately modeled, even in inhomogeneous 
systems, using the CMF method. However, one must 
keep in mind that in equilibrium, the CMF method al- 
ways underestimates the defect density on the Mott side, 
and this should be taken into account in interpreting the 
results of dynamics. 

The slow regime is dominated by long wavelength 
amplitude modes, which are described by the MF+G 
method. In this regime we find that universal power law 
scaling in the number of excitations created as a function 



of the ramp rate sets in for sufficiently slow ramps. That 
is the number of excitations (or defects for ramps that 
go deep into the Mott insulator) goes as r 2 / 3 where r is 
the ramp rate. Quantitatively, the appearance of scaling 
depends on the protocol being used. Scaling is easiest to 
see for ramps that start at the QCP and go deep into the 
Mott regime, as the scaling is not obscured by 2irh/J c 
oscillations nor the decay of the order parameter. In- 
deed, in this protocol the scaling sets in at ^(2irh/J c ) 
timescale. On the other hand, for ramps that start in SF 
phase, longer timescales of ~ 15 * (2irh/J c ) are required 
to observe 2/3 power law scaling. 



In conclusion, our calculations in equilibrium provide 
a number of results that can be useful for interpreting 
experiments, including an intuition for the behavior of 
defect density near the QCP and a simple expression for 
the defect-defect correlation function <?2(0 in the vicinity 
of the QCP. Further, understanding what our methods 
get right and what they get wrong in equilibrium gives 
us an understanding of how to apply them to dynamics. 
Our dynamical calculations provide an understanding of 
short timescales dominated by the peak in the density 
of states. Finally, we obtain the defect densities and 
timescales needed to observe universal power law scal- 
ing. 
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Appendix A: MF+G in equilibrium 



The goal of this appendix is to provide details of the modifications to the Mean Field ground state due to Gaussian 
fluctuations. We begin by giving details of the Bogoliubov transformations, and then construct expressions for the 
defect density and various correlation functions. 



1. Bogoliubov transformations 



The quadratic Hamiltonians H a ,k are solved by the Bogoliubov transformations 
where the coherence factors are 



(Al) 



Uak = 



OL<,,k 



1 



and the Hamiltonian becomes 

H " = 9 Yl ( E ° k + \ ) xtfe7<rfc 



t'ak = -sgn(/3 (Tife ) 



(A2) 



(A3) 



The ground state corresponds to the vacuum of 7 CTi fe bosons, or in terms of the b a< k operators to the squeezed state of 
Eq. f3l with 



Uk = 0. 



(A4) 



2. Particle, hole, and defect density 



To obtain an expression for the particle, hole, and defect density operators in terms of the b operators, we use the 
same procedure as we used with the Hamiltonian: we first write down the density operator in terms of the t operators 
of Eq. (16), next we replace the t operators by b operators using the transformation M from Eq. (16 1, finally we 
expand the result to second order in b a i and b^ i keeping in mind that we must replace b\ ^fro.i = 1 — " a jb a ,i ~ b^ 
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Following this procedure we obtain the operators for the particle, hole, and defect density P p (i), Ph(i), Pd{i) 



P P (i) = tl ti t +ti = \ [1 - cos(0)] - i sin(0) [&L + b ai 
+ \ cos(&)6^6 a< + - [1 + cos(<9)] - - cos 

P^i) = ti t it-, t = \ [1 - cos(0)] - J sm(0) [bl L + b al 

+ \ cos(0)?4& ai + - [1 + cos(0)] blfin + - cos ( ~ 

i , d(*) = (s?) 2 = tt, i * + ,i + tL l< t_, i 

= i [1 - cos(0)] - i sin(0) [&L + 6 m 



<?!>2 



6L6„ 



1 . 






- sin 




2 







+ cos(6')6 t Qi 6 m + - [1 + cos(0)] ft^fyi 



(A5) 



(A6) 



(A7) 



We note that the expressions we obtain for the density operators satisfy the relation that Pd(i) — P p (i) + Ph{i)- In the 
following subsection, we shall use the operators P p (i), Ph{i), Pd{i) to compute their correlation functions. However, 
in addition to the correlation functions we shall also need to compute the expectation value of the defect density 
operator P d in equilibrium in order to compare the result from MF+G method with the other methods. Taking the 
expectation value of the Pd(i) operator in equilibrium, we obtain 



(Pd) = \ [1 - cos(0)] + cos(#) Y, <h + ~ [1 + cos(0)] v . 



2 

<t>,k 



(A8) 



where the values for the coherence factors in equilibrium are obtained from expressions Eqs. (15), (28)-(31), (A2). To 
complete the discussion, we give the explicit expression for defect probability in terms of the quadratic Hamiltonian 
parameters and /3 CTi fc 



(P d (i)) = sin 



cos 



d 2 k 



a a .k 



2«/< fc -/£ft 



1 - sin 2 



d 2 k 

WT 2 



a<t>M 



(A9) 



3. Correlation functions 

In this subsection our goal is to compute the correlation functions 

f? P (l) = (PS)) 2 {gV(l) - 1) = (Pp(i)P P (i + 0) - (P P (i)){Pp(i + 0), (A10) 
ft(l) = (P^)) 2 (9 P 2 h (l) -1) = (P P (i)P h {i + l)) - (P P (i))(P h (i + l)}, (All) 
/ 2 d " d (0 - (P d (i)) 2 (gl- d (i) i) = (P d (i)P d (i + l)) - (P d (i)){P d (i + l)). (A12) 

Where, we use (P d (l)) 2 to normalize all of the gi 's so that they can be directly compared with each other. In the 
following, we shall suppress the dependence on position i since we have a homogenous system. 

We begin by obtaining expressions for the /2's up to fourth order in b operators. Although fourth order may seem 
to be overkill, we will see that the second order terms disappear on the Mott side due to sin(#) becoming zero thus 
forcing us to go to the next non-zero order. 

On the superfluid side, since sin(0) 7^ the second order term 

ft d(2) d) - ^ ou +M> (Al3) 



20 



is non-zero. Using the zero temperature expressions 

(b ak b a - k ) = (bt k bt_ k ) = u M = — ft , (A14) 



we obtain 



IT A' 



<A15) 



h jfc fct \ = ,,2 = ^ ,j (A16) 



ft d(2 \l) = ^ 1 J ^^ S (l.k)((bl k + b a ,. k )(bl_ k + b a , k )), (A17) 
sin 2 (6») /" eP/c a a>k - P a ,k , A1sn 



On the Mott side, sin(0) = and the second order term vanishes; therefore, we are forced to look to the fourth order 
term. Surprisingly, we also find that the contributions from the fourth order term are important on the superfluid 
side (see main text). We note that the only interesting fourth order term arises as a consequence of the product of 
the two second order terms in (A7) [i.e., Pjf\i) P^ii + 1) is independent of I and is therefore canceled by the regular 



part; P^\i)P^ 3 \i + I) vanishes due to sin(#) = 0]. Applying Wick's theorem, which we are allowed to do since we 
are working within the quadratic approximation to the Hamiltonian, we obtain 

ft d[4 \l) = cos 2 (#) [(bl ti bi >i+l ) (b a ,ib a , i+l ) + (bl^+i) (bl,ib a ,i+l) 



+ 



1 + cos(0) 



u <j>,i u 4>,i+i 



) (b^ib^+i) + (bl^i+ij (b\ ti b^ i+ i 



= cos 2 (9) [Fl + Gl] + - (1 + cos(#)) 2 [F 2 + G 2 ] 



where we have introduced the notation 



F„ = 



cPk j3 a ^ k cos(7 • k) 



(27T) 



G„ — 



2Jai 



d 2 k a G:k co$>{l ■ k) 



2 Jo?. 



PI 



(A19) 
(A20) 

(A21) 



By going to fourth order in fluctuations for correlations but only to second order in the Hamiltonian we miss another 
term that is fourth order in fluctuations, namely fourth order on site i and zeroth order on site j and vice versa. As 
we already explained, this term does not contribute to fa and only enters gi via a slight modification of the defect 
probability itself. 



We can find similar expressions for the particle-particle, particle-hole, and hole-hole correlation functions: 
f 2 - p{2) (l) = ft h(4) (l) = I sm 2 W [F a +G a ] + \ sin 2 (£) [F, + G,] , 
/r M2) (0 = f H-m {l) = 1 sin 2 ((?) [Fa + Ga ]-\ sin 2 (t\ [F* + Go] , 



(A22) 
(A23) 



and 



/r H4) (0 = ft m (l) = \ ™s 2 (6) [Fl + Gl] + 1 (1 + cos(0)) 2 [Fl + G 2 ] + \ cos 2 (0 [F a F, + G a G,} , (A24) 
f- h = ft m (l) = \ cos 2 (0) [Fl + Gl] + 1 (1 + cos(0)f [Fl + Gl]-\ cos 2 (0 [F a F, + G a G,} . (A25) 
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FIG. 12: Analogous to Figure [4] in the main text, but with normalized correlation functions: 

(a\ n d - d ll - 11 1 ~ ( p dWPd(0))-{Pd(0)) 2 

(a) g 2 (l-l) l = (P d( o)>2 ' 

(Vl„P-Pf7-11 i - <Pp(QPp(0))-(Pp(0)) 2 
W 52 v / 1 - <P P (0)) 2 

M n P~ h ll - 11 _ 1 = ( P p( l ) P h(0))-(Pp(0))(Ph(0)} 
(C) 92 \ L — 1 ) J- - (Pp(0)>(P h (0)> 



Appendix B: Evolution via quadratures 



An alternative view to writing down the Schrodinger equation for the wave function, Eq. (34), is to study directly 



the evolution of operators of interest using the Heisenberg equations of motion. Since the Hamiltonian is always 
quadratic, the Heisenberg equations of motion for the quadratures (b^ k b a ,k + b^ _ k b (T! - k ), (b^ ^,6^ _ fc ), and (6<r,fc&<r,fc) 
close on themselves. Taking into account the Berry phase, the Heisenberg equations of motion are 



d t(K, k ba,k 



dt{K,kK,-k) 



h 
2i 

i 



([blj.M + bl 

(-2a a , k (t)(bl k bt t _ k ) - P*, k (t)(bl ik b a , k + bl_ k b a ,. k ) 



•24, fc (*)<W,-*>' 



J (2a ff ,fe(*)(6 ff ,fefe<r,-fc) +^,fc(t)(^, fc ^,fe 



-k°<r,-k, 



+ 2A aik (t){b aik b a ._ k ). 



(Bl) 



(B2) 



(B3) 



The initial conditions for the quadratures at zero temperature can be obtained from Eqs. (A14)-(A16). A useful 



remark is that due to the quadratic nature of the effective Hamiltonian, the thermal as well as ground states of the 
system can be expressed by specifying just the three quadratures appearing in the above equations. Moreover, the 
equations of motion for the quadratures still hold for the thermal state. However, to take finite temperature into 
account we must modify the MF solution on top of which we develop the Gaussian fluctuations, i.e. 6{t) and <j>(t), as 
well as the initial conditions. 



Appendix C: Normalized correlation functions and correlation functions in dynamics 



In this appendix we supplement the main text by (1) plot ting the normalized equilibrium correlation functions, e.g. 

-d 



9~ 2 - - ft d /(Pd(i)) 2 + 1 = (Pd(i)P d (i 

(obtained using ED, CMF, and MF+G' methods) at the end of the ramp [Fig. [13 



1)) / \Pd(i)) 2 [Fig. 12 and (2) plotting the correlation function f$~ d (l = 1) 



Appendix D: Connection with Quantum Gas microscope experiments 



The absolute timescales involved in experiments depend on the details of the implementation. However, the use 
of Rb atoms is popular, and thus many setups will have similar timescales. We have used the setup of Ref [2] 
throughout the manuscript to compute typical timescales in terms of seconds. The details, which we now provide, 
that are necessary to make the connection are the dependence of the hopping matrix element J and the onsite repulsion 
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FIG. 13: Nearest neighbor correlation function f2~ d (i = 1) (immediately following the end of the ramp) as a function of the 
ramp time. The calculations were performed for four different ramp types (as indicated by row headings), using the three 
(applicable) methods (as indicated by column headings). Parameters used in calculations are identical to those of Fig. [7] 



U on the optic lattice depth V. Explicitly, we find 

J(V) = 244.5 exp(-0.209F) Hz, 
U(V) = (83.4 + 23.17 - 0.2985T/ 2 ) Hz, 

where V is measured in units of recoil energy. 
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